library(ggplot2)

boxplotToPdf <- function(file,data,stringTitle,stringXLabel,color,pdfHeight,pdfWidth){
  pdf(file)
  boxplot(as.data.frame(data),horizontal=TRUE,main=stringTitle,col=color,las=1)
  title(xlab=stringXLabel)
  dev.off()
}

extractMetrics <- function(inputFile){
  print(inputFile)
  #Read data
  inputData <- read.table(inputFile, header = TRUE)
  
  #Extract the number of robots in the chain at the last time step
  robotsInChain <- inputData[dim(inputData)[1],2]
  
  #Extract the completion time
  completionTime <- inputData[dim(inputData)[1],1]
  
  #Extract the completion time
  firstTargetTimes <- inputData[inputData[,3] > 0,1]
  
  returnValues <- c(as.numeric(robotsInChain),as.numeric(completionTime),as.numeric(firstTargetTimes[1]))
  return(returnValues)
}

padVector <- function(vector,desiredLength){
  if(length(vector) < desiredLength){
    paddingLength = desiredLength - length(vector)
    return(append(vector,rep(vector[length(vector)],paddingLength)))
  }
  return(vector)
}

extractColumn <- function(inputFile,index,maxDuration){
  #Read data
  inputData <- read.table(inputFile, header = TRUE)
  return(padVector(inputData[,index],maxDuration))
}


plotRobotInChain <- function(inChainPlot){
  #Plot max error 
  pdf("RobotInChain.pdf",family="Gentium")
  xRange <- range(0,dim(inChainPlot)[1])
  yRange <- range(0,50)
  nPlots <- dim(inChainPlot)[2]
  plotColors <- rainbow(nPlots) 
  lineType <- c(1,2,4,2,1,5) 
  plotChar <- seq(18,18+nPlots,1)
  
  
  plot(xRange,yRange, type="n", xlab="Simulation steps", ylab="Robots in chain",mar=c(7,4,1.5,2)+0.1)   
  
  #Add lines
  matlines(inChainPlot, type="l" , lwd=1.5, lty=lineType, col=plotColors, pch=plotChar)
  
  #Add a title
  #title(graphTitle)
  
  #Add a legend 
  legend(x="bottomright", colnames(inChainPlot), cex=0.8, col=plotColors, lty=lineType, title="Percentiles")
  dev.off()
  
  write.table(inChainPlot,"RobotInChainPlots.plt",sep="\t",row.names=FALSE,col.names=TRUE)
}

plotRobotOnTarget <- function(onSpotPlot){
  #Plot max error 
  pdf("RobotOnTarget.pdf",family="Gentium")
  xRange <- range(0,dim(onSpotPlot)[1])
  yRange <- range(0,7)
  nPlots <- dim(onSpotPlot)[2]
  plotColors <- rainbow(nPlots) 
  lineType <- c(1,2,4,2,1,5) 
  plotChar <- seq(18,18+nPlots,1)
  
  plot(xRange,yRange, type="n", xlab="Simulation steps", ylab="Robots on spot",mar=c(7,4,1.5,2)+0.1)   

  #Add lines
  matlines(onSpotPlot, type="l" , lwd=1.5, lty=lineType, col=plotColors, pch=plotChar)

  #Add a title
  #title(graphTitle)

  #Add a legend 
  legend(x="bottomright", colnames(onSpotPlot), cex=0.8, col=plotColors, lty=lineType, title="Percentiles")
  dev.off()

  write.table(onSpotPlot,"RobotOnSpotPlots.plt",sep="\t",row.names=FALSE,col.names=TRUE)
}


#computeQuantileByRow <- function(matrix){
#  return t(apply(matrix,1,quantile))
#}
#############################################################################
library(extrafont)

#Set the directory containing the experiment results as current working directory
setwd("Results/50Trials")
#List all the files containing information about occupation
resultsFiles <- list.files(pattern = "\\.txt$")
#Process every single file separately
print("[STATUS] - Processing results file")
metrics <- sapply(resultsFiles,extractMetrics)
#Transform processing results into dataframe
metrics <- as.matrix(metrics)
rownames(metrics) <- c("Robots in Chain","Completion Time","First Target Time")
print(metrics)
duration <- max(metrics["Completion Time",])

inChain <- lapply(resultsFiles,extractColumn,index=2,maxDuration=duration)
onTarget <- lapply(resultsFiles,extractColumn,index=3,maxDuration=duration)
inChain <- data.frame(matrix(unlist(inChain), ncol=length(inChain)))
onTarget <- data.frame(matrix(unlist(onTarget), ncol=length(onTarget)))

#print("[STATUS] - Computing quantiles & mean")
#Given the matrix containing the integrated errors for each trial as column, compute the quantiles per row.
#That corresponds to obtain the quantiles at each time step
inChainPlot <- t(apply(inChain,1,quantile))
onTargetPlot <- t(apply(onTarget,1,quantile))
#Compute row by row mean of the integrated error of every single trial
meanInChain <- apply(inChain,1,mean)
meanOnTarget <-apply(onTarget,1,mean)
#Bind toghether the elements
inChainPlot <- cbind(inChainPlot,meanInChain)
colnames(inChainPlot) <- c("0%","25%", "50%","75%","100%","Mean")
onTargetPlot <- cbind(onTargetPlot,meanOnTarget)
colnames(onTargetPlot) <- c("0%","25%", "50%","75%","100%","Mean")

print("[STATUS] - Create reporting data frame")
#Create data frame for reporting
report <- data.frame()
report["Mean","Robots in Chain"] <- mean(metrics["Robots in Chain",])
report["Mean","Completion Time"] <- mean(metrics["Completion Time",])
report["Mean","First Target Time"] <- mean(metrics["First Target Time",])
report["Standard Deviation","Robots in Chain"] <- sd(metrics["Robots in Chain",])
report["Standard Deviation","Completion Time"] <- sd(metrics["Completion Time",])
report["Standard Deviation","First Target Time"] <- sd(metrics["First Target Time",])
report["Median","Robots in Chain"] <- median(metrics["Robots in Chain",])
report["Median","Completion Time"] <- median(metrics["Completion Time",])
report["Median","First Target Time"] <- median(metrics["First Target Time",])
report["Min","Robots in Chain"] <- min(metrics["Robots in Chain",])
report["Min","Completion Time"] <- min(metrics["Completion Time",])
report["Min","First Target Time"] <- min(metrics["First Target Time",])
report["Max","Robots in Chain"] <- max(metrics["Robots in Chain",])
report["Max","Completion Time"] <- max(metrics["Completion Time",])
report["Max","First Target Time"] <- max(metrics["First Target Time",])
#Write reporting dataframe
write.table(report,"Results.stat",sep="\t",row.names=TRUE,col.names=TRUE)
#write(paste("Correlation (Pearson)",cor(x=metrics["Robots in Chain",],y=metrics["Completion Time",], use="all.obs", method="pearson"),sep="\t"),"Result.stat",append = TRUE)
#write(paste("Correlation (Spearman)",cor(x=metrics["Robots in Chain",],y=metrics["Completion Time",], use="all.obs", method="spearman"),sep="\t"),"Result.stat",append = TRUE)
#write(paste("Correlation (Kendall)",cor(x=metrics["Robots in Chain",],y=metrics["Completion Time",], use="all.obs", method="kendall"),sep="\t"),"Result.stat",append = TRUE)

#Plot box plots and empirical cumulative density functions
#boxplotToPdf("BoxplotTime.pdf",metrics["Completion Time",],paste("Boxplot of completion times across",length(resultsFiles),"trials"),"Simulation steps","blue")
#boxplotToPdf("BoxplotRobots.pdf",metrics["Robots in Chain",],paste("Boxplot of number of robots in chain across",length(resultsFiles),"trials"),"Number of robots","orange")

# Add boxplots to a scatterplot
# Scatterplot
pdf("ResultsDistribution.pdf", family="Gentium")
par(fig=c(0,0.8,0,0.8), new=TRUE, bty="l")
plot(x=metrics["Robots in Chain",],y=metrics["Completion Time",],col=rgb(0,100,0,80,maxColorValue=255), pch=19 , xlab="Robots in Chain",
     ylab="Completion Time",cex=0.8,cex.lab=1.5)
# Fit lines
abline(lm(metrics["Completion Time",]~metrics["Robots in Chain",]), col="red") # regression line (y~x) 
#lines(lowess(metrics["Robots in Chain",],metrics["Completion Time",]), col="blue") # lowess line (x,y)
#legend(x="bottomright", c(), cex=0.8, col=plotColors, lty=lineType, title="Lines")
# Boxplot X-Axis
par(fig=c(0,0.8,0.55,1), new=TRUE)
boxplot(metrics["Robots in Chain",], horizontal=TRUE, axes=FALSE)
# Boxplot Y-Axis
par(fig=c(0.65,1,0,0.8), new=TRUE)
boxplot(metrics["Completion Time",], axes=FALSE)
# Title
#mtext(expression(bold("Scatterplot of robots in chain vs completion time")), side=3, outer=TRUE, line=-3,cex=1.0)
dev.off()

#Plot empirical CDF of completion times
pdf("Distributions.pdf",family="Gentium")
par(mfrow=c(2,2))
hist(metrics["Completion Time",],col="blue",xlab="Simulation steps", main=paste("Histogram of completion times\nacross",length(resultsFiles),"trials"))
plot(ecdf(metrics["Completion Time",]),col="blue",xlab="Simulation steps",main=paste("Empirical cdf of completion times\nacross",length(resultsFiles),"trials"))
hist(metrics["Robots in Chain",],col="orange",xlab="Robots",main=paste("Histogram of number of robots in chain\nacross",length(resultsFiles),"trials"))
plot(ecdf(metrics["Robots in Chain",]),col="orange",xlab="Robots",main=paste("Empirical cdf number of robots in chain\nacross",length(resultsFiles),"trials"))
dev.off()

pdf("Panel1.pdf",family="Gentium")
hist(metrics["Completion Time",],col="blue",xlab="Simulation steps", mar=c(7,4,1.5,2)+0.1, cex.lab=1.5, main="")
dev.off()

pdf("Panel2.pdf",family="Gentium")
plot(ecdf(metrics["Completion Time",]),col="blue",xlab="Simulation steps",mar=c(7,4,1.5,2)+0.1, cex.lab=1.5, main="")
dev.off()

pdf("Panel3.pdf",family="Gentium")
hist(metrics["Robots in Chain",],col="orange",xlab="Robots",mar=c(7,4,1.5,2)+0.1, cex.lab=1.5, main="")
dev.off()

pdf("Panel4.pdf",family="Gentium")
plot(ecdf(metrics["Robots in Chain",]),col="orange",xlab="Robots",mar=c(7,4,1.5,2)+0.1, cex.lab=1.5, main="")
dev.off()



#plotRobotInChain(inChainPlot)
#plotRobotOnTarget(onTargetPlot)

print("Correlation - Pearson")
print(cor(x=metrics["Robots in Chain",],y=metrics["Completion Time",],method="pearson"))
print(cor.test(x=metrics["Robots in Chain",],y=metrics["Completion Time",],method="pearson"))

print("Correlation - Kendall")
print(cor(x=metrics["Robots in Chain",],y=metrics["Completion Time",],method="kendall"))
print(cor.test(x=metrics["Robots in Chain",],y=metrics["Completion Time",],method="kendall"))

print("Correlation - Spearman")
print(cor(x=metrics["Robots in Chain",],y=metrics["Completion Time",],method="spearman"))
print(cor.test(x=metrics["Robots in Chain",],y=metrics["Completion Time",],method="spearman"))
#Output analysis result
#sprintf("Time required to get 25%% coverage: %.0f",firstQuartileOccupation[1])
#sprintf("Time required to get 50%% coverage: %.0f",secondQuartileOccupation[1])
#sprintf("Time required to get 75%% coverage: %.0f",thirdQuartileOccupation[1])

rm(metrics)
rm(report)
rm(resultsFiles)

#createTexGraph("t","$\\max_i(e_i(t))$",5,5,inputData[,1],maxErrorIntegrated[,1])
